1 Setup and data wrangling

reading files and adding/checking LTER names

sisyn <- read_csv("20201111_masterdata_RAW.csv") 

wrtds_annual <- read_csv("WRTDS_AnnualResults_AllSites_052621.csv") #%>% 
  #rename(site=Site)

wrtds_trends <- read_csv("Si_EGRETCi_Trends_AllSites.csv") %>% 
  rename(site=Site)



sisyn_sites <- sisyn %>%  select(site, LTER) %>% unique()

wrtds_annual <- left_join(wrtds_annual, sisyn_sites)

wrtds_trends <- left_join(wrtds_trends, sisyn_sites) %>% 
  mutate(LTER_site=paste(LTER, site, sep = "------"))

x <- wrtds_annual %>% filter(is.na(LTER))
unique(x$site)
## [1] "ws1" "ws2" "ws3" "ws4" "ws5" "ws6" "ws7" "ws8" "ws9"
wrtds_annual <- wrtds_annual %>% 
  mutate(LTER=ifelse(is.na(LTER), "HBR", LTER))

wrtds_trends <- wrtds_trends %>% 
  mutate(LTER=ifelse(is.na(LTER), "HBR", LTER))

z <- wrtds_annual %>% count(LTER,site)
z
## # A tibble: 52 x 3
##    LTER  site              n
##    <chr> <chr>         <int>
##  1 AND   GSMACK           24
##  2 AND   GSWS02           37
##  3 AND   GSWS06           17
##  4 AND   GSWS07           18
##  5 AND   GSWS08           37
##  6 AND   GSWS09           37
##  7 AND   GSWS10           37
##  8 ARC   Imnavait Weir    17
##  9 HBR   ws1              48
## 10 HBR   ws2              54
## # ... with 42 more rows

pivoting the annual discharge/flux/conc/FNconc data and then calculating correlations and regressions across time

wrtds_long <- wrtds_annual %>%
  pivot_longer(cols = c(Discharge_cms:FNFlux_106_kg_y),
               names_to = "variable",
               values_to = "value")

#correlations:
wrtds_corr <- wrtds_long %>% 
  group_by(site, LTER, variable) %>% 
  summarize(ri=cor(Year, value, use = "pairwise.complete.obs"), ni=n())


#regressions:
wrtds_slopes <-  wrtds_long %>%
  filter(!is.na(value)) %>% 
  group_by(site, LTER, variable) %>%
  do(model = tidy(lm(value ~ Year, data = .))) %>% # only works if you include the broom::tidy in here
  unnest(model)

2 Correlations/slopes from annual estimates

ggplot(wrtds_corr, aes(LTER, ri)) + 
    stat_summary(fun = mean, geom = "point") + 
    stat_summary(fun.data = mean_cl_normal, geom = "pointrange", fun.args = list(mult = 1))+
  geom_hline(yintercept = 0)+
  xlab("LTER")+
  ylab("correlation of DSi ~ year")+
  stat_summary(fun.data = give.n, geom = "text", size=3, fontface="italic", position = position_nudge(x = 0.4, y=.1))+
  coord_flip()+
  scale_x_discrete(limits = rev)+
  facet_wrap(~variable)

ggplot(wrtds_corr, aes(site, ri)) + 
    stat_summary(fun = mean, geom = "point") + 
    stat_summary(fun.data = mean_cl_normal, geom = "pointrange", fun.args = list(mult = 1))+
  geom_hline(yintercept = 0)+
  xlab("LTER")+
  ylab("correlation of DSi ~ year")+
  stat_summary(fun.data = give.n, geom = "text", size=3, fontface="italic", position = position_nudge(x = 0.4, y=.1))+
  coord_flip()+
  scale_x_discrete(limits = rev)+
  facet_wrap(~variable)

discharge is super correlated to flux but not conc across sites and years? does this make sense?

pairs.panels(wrtds_annual)

models using correlation data

ES_wrtds <- escalc(measure="ZCOR", ri=ri, ni=ni, data=wrtds_corr)

mod1<-rma(yi=yi,vi=vi,data=filter(ES_wrtds, variable=='Discharge_cms'))
summary(mod1)
## 
## Random-Effects Model (k = 52; tau^2 estimator: REML)
## 
##   logLik  deviance       AIC       BIC      AICc 
##   1.3614   -2.7228    1.2772    5.1409    1.5272   
## 
## tau^2 (estimated amount of total heterogeneity): 0.0089 (SE = 0.0098)
## tau (square root of estimated tau^2 value):      0.0944
## I^2 (total heterogeneity / total variability):   17.62%
## H^2 (total variability / sampling variability):  1.21
## 
## Test for Heterogeneity:
## Q(df = 51) = 63.2050, p-val = 0.1173
## 
## Model Results:
## 
## estimate      se    zval    pval   ci.lb   ci.ub 
##   0.1626  0.0315  5.1603  <.0001  0.1009  0.2244  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod2<-rma(yi=yi,vi=vi,data=filter(ES_wrtds, variable=='Conc_mgL'))
summary(mod2)
## 
## Random-Effects Model (k = 52; tau^2 estimator: REML)
## 
##   logLik  deviance       AIC       BIC      AICc 
## -49.4216   98.8432  102.8432  106.7069  103.0932   
## 
## tau^2 (estimated amount of total heterogeneity): 0.3623 (SE = 0.0809)
## tau (square root of estimated tau^2 value):      0.6019
## I^2 (total heterogeneity / total variability):   89.68%
## H^2 (total variability / sampling variability):  9.69
## 
## Test for Heterogeneity:
## Q(df = 51) = 528.5065, p-val < .0001
## 
## Model Results:
## 
## estimate      se    zval    pval    ci.lb   ci.ub 
##   0.0619  0.0887  0.6985  0.4849  -0.1119  0.2357    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod3<-rma(yi=yi,vi=vi,data=filter(ES_wrtds, variable=='FNConc_mgL'))
summary(mod3)
## 
## Random-Effects Model (k = 52; tau^2 estimator: REML)
## 
##   logLik  deviance       AIC       BIC      AICc 
## -79.7092  159.4184  163.4184  167.2821  163.6684   
## 
## tau^2 (estimated amount of total heterogeneity): 1.2860 (SE = 0.2639)
## tau (square root of estimated tau^2 value):      1.1340
## I^2 (total heterogeneity / total variability):   96.86%
## H^2 (total variability / sampling variability):  31.83
## 
## Test for Heterogeneity:
## Q(df = 51) = 1529.9495, p-val < .0001
## 
## Model Results:
## 
## estimate      se     zval    pval    ci.lb   ci.ub 
##  -0.0859  0.1601  -0.5363  0.5917  -0.3996  0.2279    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod3.1<-rma(yi=yi,vi=vi,data=filter(ES_wrtds, variable=='FNConc_mgL'),
            mods = ~LTER)
summary(mod3.1)
## 
## Mixed-Effects Model (k = 52; tau^2 estimator: REML)
## 
##   logLik  deviance       AIC       BIC      AICc 
## -62.3706  124.7411  146.7411  165.8555  155.5411   
## 
## tau^2 (estimated amount of residual heterogeneity):     1.0938 (SE = 0.2487)
## tau (square root of estimated tau^2 value):             1.0459
## I^2 (residual heterogeneity / unaccounted variability): 96.37%
## H^2 (unaccounted variability / sampling variability):   27.52
## R^2 (amount of heterogeneity accounted for):            14.94%
## 
## Test for Residual Heterogeneity:
## QE(df = 42) = 1096.5069, p-val < .0001
## 
## Test of Moderators (coefficients 2:10):
## QM(df = 9) = 17.5632, p-val = 0.0406
## 
## Model Results:
## 
##                        estimate      se     zval    pval    ci.lb    ci.ub 
## intrcpt                  0.0112  0.4030   0.0277  0.9779  -0.7787   0.8010    
## LTERARC                  1.3077  1.1522   1.1349  0.2564  -0.9507   3.5660    
## LTERHBR                  0.3572  0.5361   0.6663  0.5052  -0.6935   1.4079    
## LTERKRR(Julian)          0.0673  0.6713   0.1003  0.9201  -1.2483   1.3830    
## LTERLMP(Wymore)          0.8226  1.1546   0.7124  0.4762  -1.4405   3.0856    
## LTERLUQ                  0.4282  0.7355   0.5822  0.5604  -1.0133   1.8698    
## LTERMCM                 -1.3196  0.5525  -2.3885  0.0169  -2.4024  -0.2368  * 
## LTERNWT                  0.7491  0.7376   1.0155  0.3098  -0.6966   2.1948    
## LTERSagehen(Sullivan)   -0.4160  1.1385  -0.3654  0.7148  -2.6475   1.8154    
## LTERUMR(Jankowski)      -0.2156  0.4886  -0.4413  0.6590  -1.1733   0.7420    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

2.1 meta-regression models

wrtds_slopes <- left_join(wrtds_slopes, wrtds_corr) %>% 
  mutate(mi=1) %>% 
  filter(term=='Year')

wrtds_slopes <- wrtds_slopes %>% 
  arrange(LTER)

wrtds_slopes_ES_Si <- escalc(measure='ZPCOR',
                          ti=statistic,
                          ni=ni,
                          mi=mi,
                          data=filter(wrtds_slopes, variable=='Conc_mgL'))

wrtds_slopes_ES_discharge <- escalc(measure='ZPCOR',
                          ti=statistic,
                          ni=ni,
                          mi=mi,
                          data=filter(wrtds_slopes, variable=='Discharge_cms'))

wrtds_slopes_ES_flux <- escalc(measure='ZPCOR',
                          ti=statistic,
                          ni=ni,
                          mi=mi,
                          data=filter(wrtds_slopes, variable=='Flux_106_kg_y'))


mod_Si1<-rma.mv(yi,vi,
          data=filter(wrtds_slopes_ES_Si),
          random = ~1|site)
summary(mod_Si1)
## 
## Multivariate Meta-Analysis Model (k = 52; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
## -49.4102   98.8205  102.8205  106.6841  103.0705   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed  factor 
## sigma^2    0.3641  0.6034     52     no    site 
## 
## Test for Heterogeneity:
## Q(df = 51) = 549.1464, p-val < .0001
## 
## Model Results:
## 
## estimate      se    zval    pval    ci.lb   ci.ub 
##   0.0621  0.0886  0.7005  0.4836  -0.1116  0.2358    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
forest(mod_Si1)

mod_Si2<-rma.mv(yi,vi,
          data=filter(wrtds_slopes_ES_Si))
summary(mod_Si2)
## 
## Multivariate Meta-Analysis Model (k = 52; method: REML)
## 
##    logLik   Deviance        AIC        BIC       AICc 
## -240.9324   481.8647   483.8647   485.7965   483.9464   
## 
## Variance Components: none
## 
## Test for Heterogeneity:
## Q(df = 51) = 549.1464, p-val < .0001
## 
## Model Results:
## 
## estimate      se    zval    pval    ci.lb   ci.ub 
##   0.0523  0.0277  1.8867  0.0592  -0.0020  0.1066  . 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
forest(mod_Si2)

mod_discharge1<-rma.mv(yi,vi,
          data=filter(wrtds_slopes_ES_discharge),
          random = ~1|site)
summary(mod_discharge1)
## 
## Multivariate Meta-Analysis Model (k = 52; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
##   1.3103   -2.6207    1.3793    5.2430    1.6293   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed  factor 
## sigma^2    0.0109  0.1043     52     no    site 
## 
## Test for Heterogeneity:
## Q(df = 51) = 66.1827, p-val = 0.0749
## 
## Model Results:
## 
## estimate      se    zval    pval   ci.lb   ci.ub 
##   0.1628  0.0316  5.1462  <.0001  0.1008  0.2248  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
forest(mod_discharge1)

mod_flux1<-rma.mv(yi,vi,
          data=filter(wrtds_slopes_ES_flux),
          random = ~1|site)
summary(mod_flux1)
## 
## Multivariate Meta-Analysis Model (k = 52; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
##  -5.8999   11.7998   15.7998   19.6635   16.0498   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed  factor 
## sigma^2    0.0335  0.1831     52     no    site 
## 
## Test for Heterogeneity:
## Q(df = 51) = 93.2029, p-val = 0.0003
## 
## Model Results:
## 
## estimate      se    zval    pval   ci.lb   ci.ub 
##   0.2005  0.0382  5.2491  <.0001  0.1256  0.2754  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
forest(mod_flux1)

viz_forest(x = wrtds_slopes_ES_Si[, c("yi", "vi")], 
           study_labels = wrtds_slopes_ES_Si[, "site"],
           xlab = "Slope of Si by year (Z-transformed)",
           variant = "thick",
           group=wrtds_slopes_ES_Si[, "LTER"],
           method = "REML",
           annotate_CI = TRUE,
           summary_label = c("Summary(LTER = AND)",
                  "Summary(LTER = ARC)",
                  "Summary(LTER = KRR(Julian))",
                  "Summary(LTER = HBR)",
                  "Summary(LTER = LMP(Wymore))",
                  "Summary(LTER = LUQ)",
                  "Summary (LTER = MCM)",
                  "Summary (LTER = NWT)",
                   "Summary(LTER = Sagehen(Sullivan))",
                   "Summary(LTER = UMR(Jankowski))"))

2.2 forest plots for modeled DSi conc by LTER

a <- viz_forest(x = filter(wrtds_slopes_ES_Si, LTER=="AND")[, c("yi", "vi")], 
           study_labels = filter(wrtds_slopes_ES_Si, LTER=="AND")[, "site"],
           xlab = "ANDREWS: Slope of DSi by year (Z-transf.)",
           variant = "rain")
a     

b <- viz_forest(x = filter(wrtds_slopes_ES_Si, LTER=="ARC")[, c("yi", "vi")], 
           study_labels = filter(wrtds_slopes_ES_Si, LTER=="ARC")[, "site"],
           xlab = "ARCTIC: Slope of DSi by year (Z-transf.)",
           variant = "rain")
b    

c <- viz_forest(x = filter(wrtds_slopes_ES_Si, LTER=="HBR")[, c("yi", "vi")], 
           study_labels = filter(wrtds_slopes_ES_Si, LTER=="HBR")[, "site"],
           xlab = "HUBBARD BROOK: Slope of DSi by year (Z-transf.)",
           variant = "rain")
c

d <- viz_forest(x = filter(wrtds_slopes_ES_Si, LTER=="KRR(Julian)")[, c("yi", "vi")], 
           study_labels = filter(wrtds_slopes_ES_Si, LTER=="KRR(Julian)")[, "site"],
           xlab = "KISSIMMEE: Slope of DSi by year (Z-transf.)",
           variant = "rain")
d

e <- viz_forest(x = filter(wrtds_slopes_ES_Si, LTER=="LMP(Wymore)")[, c("yi", "vi")], 
           study_labels = filter(wrtds_slopes_ES_Si, LTER=="LMP(Wymore)")[, "site"],
           xlab = "LAMPREY: Slope of DSi by year (Z-transf.)",
           variant = "rain")
e

f <- viz_forest(x = filter(wrtds_slopes_ES_Si, LTER=="LUQ")[, c("yi", "vi")], 
           study_labels = filter(wrtds_slopes_ES_Si, LTER=="LUQ")[, "site"],
           xlab = "LUQUILLO: Slope of DSi by year (Z-transf.)",
           variant = "rain")
f

g <- viz_forest(x = filter(wrtds_slopes_ES_Si, LTER=="MCM")[, c("yi", "vi")], 
           study_labels = filter(wrtds_slopes_ES_Si, LTER=="MCM")[, "site"],
           xlab = "MCMURDO: Slope of DSi by year (Z-transf.)",
           variant = "rain")
g

h <- viz_forest(x = filter(wrtds_slopes_ES_Si, LTER=="NWT")[, c("yi", "vi")], 
           study_labels = filter(wrtds_slopes_ES_Si, LTER=="NWT")[, "site"],
           xlab = "NIWOT: Slope of DSi by year (Z-transf.)",
           variant = "rain")
h

i <- viz_forest(x = filter(wrtds_slopes_ES_Si, LTER=="Sagehen(Sullivan)")[, c("yi", "vi")], 
           study_labels = filter(wrtds_slopes_ES_Si, LTER=="Sagehen(Sullivan)")[, "site"],
           xlab = "Sagehen: Slope of DSi by year (Z-transf.)",
           variant = "rain")
i

j <- viz_forest(x = filter(wrtds_slopes_ES_Si, LTER=="UMR(Jankowski)")[, c("yi", "vi")], 
           study_labels = filter(wrtds_slopes_ES_Si, LTER=="UMR(Jankowski)")[, "site"],
           xlab = "Sagehen: Slope of DSi by year (Z-transf.)",
           variant = "rain")
j

ggarrange(a,b,c,d,e,f,g,h,i,j, ncol=5, nrow=2)

viz_forest(x = wrtds_slopes_ES_Si[, c("yi", "vi")], 
           xlab = "Slope of DSi by year (Z-transformed)",
           group=wrtds_slopes_ES_Si[, "LTER"],
           variant = "rain",
           study_labels = wrtds_slopes_ES_Si[, "site"],
           summary_label = c("Summary(LTER = AND)",
                  "Summary(LTER = ARC)",
                  "Summary(LTER = KRR(Julian))",
                  "Summary(LTER = HBR)",
                  "Summary(LTER = LMP(Wymore))",
                  "Summary(LTER = LUQ)",
                  "Summary (LTER = MCM)",
                  "Summary (LTER = NWT)",
                   "Summary(LTER = Sagehen(Sullivan))",
                   "Summary(LTER = UMR(Jankowski))"))

2.3 forest plots for discharge

a <- viz_forest(x = filter(wrtds_slopes_ES_discharge, LTER=="AND")[, c("yi", "vi")], 
           study_labels = filter(wrtds_slopes_ES_discharge, LTER=="AND")[, "site"],
           xlab = "ANDREWS: Slope of discharge by year (Z-transf.)",
           variant = "rain")
a     

b <- viz_forest(x = filter(wrtds_slopes_ES_discharge, LTER=="ARC")[, c("yi", "vi")], 
           study_labels = filter(wrtds_slopes_ES_discharge, LTER=="ARC")[, "site"],
           xlab = "ARCTIC: Slope of discharge by year (Z-transf.)",
           variant = "rain")
b    

c <- viz_forest(x = filter(wrtds_slopes_ES_discharge, LTER=="HBR")[, c("yi", "vi")], 
           study_labels = filter(wrtds_slopes_ES_discharge, LTER=="HBR")[, "site"],
           xlab = "HUBBARD BROOK: Slope of discharge by year (Z-transf.)",
           variant = "rain")
c

d <- viz_forest(x = filter(wrtds_slopes_ES_discharge, LTER=="KRR(Julian)")[, c("yi", "vi")], 
           study_labels = filter(wrtds_slopes_ES_discharge, LTER=="KRR(Julian)")[, "site"],
           xlab = "KISSIMMEE: Slope of discharge by year (Z-transf.)",
           variant = "rain")
d

e <- viz_forest(x = filter(wrtds_slopes_ES_discharge, LTER=="LMP(Wymore)")[, c("yi", "vi")], 
           study_labels = filter(wrtds_slopes_ES_discharge, LTER=="LMP(Wymore)")[, "site"],
           xlab = "LAMPREY: Slope of discharge by year (Z-transf.)",
           variant = "rain")
e

f <- viz_forest(x = filter(wrtds_slopes_ES_discharge, LTER=="LUQ")[, c("yi", "vi")], 
           study_labels = filter(wrtds_slopes_ES_discharge, LTER=="LUQ")[, "site"],
           xlab = "LUQUILLO: Slope of discharge by year (Z-transf.)",
           variant = "rain")
f

g <- viz_forest(x = filter(wrtds_slopes_ES_discharge, LTER=="MCM")[, c("yi", "vi")], 
           study_labels = filter(wrtds_slopes_ES_discharge, LTER=="MCM")[, "site"],
           xlab = "MCMURDO: Slope of discharge by year (Z-transf.)",
           variant = "rain")
g

h <- viz_forest(x = filter(wrtds_slopes_ES_discharge, LTER=="NWT")[, c("yi", "vi")], 
           study_labels = filter(wrtds_slopes_ES_discharge, LTER=="NWT")[, "site"],
           xlab = "NIWOT: Slope of discharge by year (Z-transf.)",
           variant = "rain")
h

i <- viz_forest(x = filter(wrtds_slopes_ES_discharge, LTER=="Sagehen(Sullivan)")[, c("yi", "vi")], 
           study_labels = filter(wrtds_slopes_ES_discharge, LTER=="Sagehen(Sullivan)")[, "site"],
           xlab = "Sagehen: Slope of discharge by year (Z-transf.)",
           variant = "rain")
i

j <- viz_forest(x = filter(wrtds_slopes_ES_discharge, LTER=="UMR(Jankowski)")[, c("yi", "vi")], 
           study_labels = filter(wrtds_slopes_ES_discharge, LTER=="UMR(Jankowski)")[, "site"],
           xlab = "Sagehen: Slope of discharge by year (Z-transf.)",
           variant = "rain")
j

ggarrange(a,b,c,d,e,f,g,h,i,j, ncol=5, nrow=2)

viz_forest(x = wrtds_slopes_ES_discharge[, c("yi", "vi")], 
           xlab = "Slope of discharge by year (Z-transformed)",
           group=wrtds_slopes_ES_discharge[, "LTER"],
           study_labels = wrtds_slopes_ES_discharge[, "site"],
           variant = "rain",
           summary_label = c("Summary(LTER = AND)",
                  "Summary(LTER = ARC)",
                  "Summary(LTER = KRR(Julian))",
                  "Summary(LTER = HBR)",
                  "Summary(LTER = LMP(Wymore))",
                  "Summary(LTER = LUQ)",
                  "Summary (LTER = MCM)",
                  "Summary (LTER = NWT)",
                   "Summary(LTER = Sagehen(Sullivan))",
                   "Summary(LTER = UMR(Jankowski))"))

2.4 forest plots for flux

a <- viz_forest(x = filter(wrtds_slopes_ES_flux, LTER=="AND")[, c("yi", "vi")], 
           study_labels = filter(wrtds_slopes_ES_flux, LTER=="AND")[, "site"],
           xlab = "ANDREWS: slope of DSi flux by year (Z-transf.)",
           variant = "rain")
a     

b <- viz_forest(x = filter(wrtds_slopes_ES_flux, LTER=="ARC")[, c("yi", "vi")], 
           study_labels = filter(wrtds_slopes_ES_flux, LTER=="ARC")[, "site"],
           xlab = "ARCTIC: slope of DSi flux by year (Z-transf.)",
           variant = "rain")
b    

c <- viz_forest(x = filter(wrtds_slopes_ES_flux, LTER=="HBR")[, c("yi", "vi")], 
           study_labels = filter(wrtds_slopes_ES_flux, LTER=="HBR")[, "site"],
           xlab = "HUBBARD BROOK: slope of DSi flux by year (Z-transf.)",
           variant = "rain")
c

d <- viz_forest(x = filter(wrtds_slopes_ES_flux, LTER=="KRR(Julian)")[, c("yi", "vi")], 
           study_labels = filter(wrtds_slopes_ES_flux, LTER=="KRR(Julian)")[, "site"],
           xlab = "KISSIMMEE: slope of DSi flux by year (Z-transf.)",
           variant = "rain")
d

e <- viz_forest(x = filter(wrtds_slopes_ES_flux, LTER=="LMP(Wymore)")[, c("yi", "vi")], 
           study_labels = filter(wrtds_slopes_ES_flux, LTER=="LMP(Wymore)")[, "site"],
           xlab = "LAMPREY: slope of DSi flux by year (Z-transf.)",
           variant = "rain")
e

f <- viz_forest(x = filter(wrtds_slopes_ES_flux, LTER=="LUQ")[, c("yi", "vi")], 
           study_labels = filter(wrtds_slopes_ES_flux, LTER=="LUQ")[, "site"],
           xlab = "LUQUILLO: slope of DSi flux by year (Z-transf.)",
           variant = "rain")
f

g <- viz_forest(x = filter(wrtds_slopes_ES_flux, LTER=="MCM")[, c("yi", "vi")], 
           study_labels = filter(wrtds_slopes_ES_flux, LTER=="MCM")[, "site"],
           xlab = "MCMURDO: slope of DSi flux by year (Z-transf.)",
           variant = "rain")
g

h <- viz_forest(x = filter(wrtds_slopes_ES_flux, LTER=="NWT")[, c("yi", "vi")], 
           study_labels = filter(wrtds_slopes_ES_flux, LTER=="NWT")[, "site"],
           xlab = "NIWOT: slope of DSi flux by year (Z-transf.)",
           variant = "rain")
h

i <- viz_forest(x = filter(wrtds_slopes_ES_flux, LTER=="Sagehen(Sullivan)")[, c("yi", "vi")], 
           study_labels = filter(wrtds_slopes_ES_flux, LTER=="Sagehen(Sullivan)")[, "site"],
           xlab = "Sagehen: slope of DSi flux by year (Z-transf.)",
           variant = "rain")
i

j <- viz_forest(x = filter(wrtds_slopes_ES_flux, LTER=="UMR(Jankowski)")[, c("yi", "vi")], 
           study_labels = filter(wrtds_slopes_ES_flux, LTER=="UMR(Jankowski)")[, "site"],
           xlab = "Sagehen: slope of DSi flux by year",
           variant = "rain")
j

ggarrange(a,b,c,d,e,f,g,h,i,j, ncol=5, nrow=2)

viz_forest(x = wrtds_slopes_ES_flux[, c("yi", "vi")], 
           xlab = "slope of DSi flux by year (Z-transformed)",
           group=wrtds_slopes_ES_flux[, "LTER"],
           study_labels = wrtds_slopes_ES_flux[, "site"],
           variant = "rain",
           summary_label = c("Summary(LTER = AND)",
                  "Summary(LTER = ARC)",
                  "Summary(LTER = KRR(Julian))",
                  "Summary(LTER = HBR)",
                  "Summary(LTER = LMP(Wymore))",
                  "Summary(LTER = LUQ)",
                  "Summary (LTER = MCM)",
                  "Summary (LTER = NWT)",
                   "Summary(LTER = Sagehen(Sullivan))",
                   "Summary(LTER = UMR(Jankowski))"))

# fit <- nlme(data=wrtds, FNConc_mgL ~ Year, random = ~LTER)
# 
# fit <- lmer(data=wrtds2, FNConc_mgL ~ Year +1|LTER)
# summary(fit)

3 Using WRTDS trend estimates

hmm some have 95% CIs in the 90 million-ish range so let’s scrap those for now

ggplot(wrtds_trends, aes(LTER_site, estC)) +
  geom_point()+
  geom_linerange(aes(ymin = lowC95, ymax = upC95))+
  coord_flip()+
  geom_hline(yintercept = 0)+
  scale_x_discrete(limits = rev)

wrtds_trends %>% 
  arrange(lowC95)
## # A tibble: 52 x 33
##    LTER  site  rejectC  pValC   estC     lowC    upC   lowC50  upC50  lowC95
##    <chr> <chr> <lgl>    <dbl>  <dbl>    <dbl>  <dbl>    <dbl>  <dbl>   <dbl>
##  1 UMR   M078~ FALSE   0.146  -0.446 -9.46e+7  0.236 -8.71e-1 -0.449 -9.53e7
##  2 UMR_~ CH00~ FALSE   0.503   0.305 -4.94e-1  1.15  -2.63e-3  0.720 -1.47e6
##  3 UMR_~ CU11~ FALSE   0.271  -0.771 -2.00e+0  0.743 -1.23e+0 -0.358 -2.36e5
##  4 UMR_~ LM00~ FALSE   0.584   0.146 -7.68e-1  0.906 -1.24e-1  0.536 -1.67e5
##  5 LUQ   RI    FALSE   0.886   1.33  -2.28e+4  5.66  -1.69e+4  2.16  -2.55e4
##  6 Sage~ Sage~ TRUE    0.0198 -1.41  -2.14e+4 -0.720 -2.09e+4 -1.21  -2.18e4
##  7 NWT   MART~ TRUE    0.0944  0.263  3.73e-2  0.414  1.65e-1  0.325 -1.77e4
##  8 ARC   Imna~ FALSE   0.795   0.553 -6.63e+0 13.1   -5.95e-1  0.718 -3.69e1
##  9 NWT   SADD~ FALSE   0.374   5.34  -5.51e+0 10.2    1.23e+0  6.97  -1.75e1
## 10 LUQ   MPR   FALSE   0.834  -0.246 -4.98e+0  3.25  -2.94e+0  0.786 -6.74e0
## # ... with 42 more rows, and 23 more variables: upC95 <dbl>, likeCUp <dbl>,
## #   likeCDown <dbl>, rejectF <lgl>, pValF <dbl>, estF <dbl>, lowF <dbl>,
## #   upF <dbl>, lowF50 <dbl>, upF50 <dbl>, lowF95 <dbl>, upF95 <dbl>,
## #   likeFUp <dbl>, likeFDown <dbl>, baseConc <dbl>, baseFlux <dbl>,
## #   iBoot <dbl>, startSeed <dbl>, nBootGood <dbl>, Year.1 <dbl>, Year.2 <dbl>,
## #   duration <dbl>, LTER_site <chr>
wrtds_trends %>% 
  arrange(desc(upC95))
## # A tibble: 52 x 33
##    LTER  site  rejectC  pValC   estC     lowC   upC   lowC50  upC50   lowC95
##    <chr> <chr> <lgl>    <dbl>  <dbl>    <dbl> <dbl>    <dbl>  <dbl>    <dbl>
##  1 ARC   Imna~ FALSE   0.795   0.553 -6.63e+0 13.1  -5.95e-1  0.718 -3.69e+1
##  2 NWT   SADD~ FALSE   0.374   5.34  -5.51e+0 10.2   1.23e+0  6.97  -1.75e+1
##  3 LUQ   RI    FALSE   0.886   1.33  -2.28e+4  5.66 -1.69e+4  2.16  -2.55e+4
##  4 UMR_~ WP02~ FALSE   0.433  -1.52  -2.10e+0  1.44 -1.63e+0 -0.157 -2.37e+0
##  5 LUQ   MPR   FALSE   0.834  -0.246 -4.98e+0  3.25 -2.94e+0  0.786 -6.74e+0
##  6 MCM   Ande~ FALSE   0.339  -0.309 -6.11e-1  2.30 -4.52e-1 -0.120 -7.45e-1
##  7 UMR_~ BK01~ FALSE   0.138   0.482 -9.85e-2  1.69  2.74e-1  1.04  -3.44e-1
##  8 KRR   S65E  FALSE   0.602   0.280 -2.36e+0  1.52 -2.62e-1  1.01  -3.18e+0
##  9 UMR_~ CN00~ FALSE   0.504   0.223 -4.01e-1  1.57 -1.16e-2  0.970 -6.30e-1
## 10 LUQ   QS    TRUE    0.0392  0.990  4.70e-1  1.64  7.72e-1  1.26   2.13e-1
## # ... with 42 more rows, and 23 more variables: upC95 <dbl>, likeCUp <dbl>,
## #   likeCDown <dbl>, rejectF <lgl>, pValF <dbl>, estF <dbl>, lowF <dbl>,
## #   upF <dbl>, lowF50 <dbl>, upF50 <dbl>, lowF95 <dbl>, upF95 <dbl>,
## #   likeFUp <dbl>, likeFDown <dbl>, baseConc <dbl>, baseFlux <dbl>,
## #   iBoot <dbl>, startSeed <dbl>, nBootGood <dbl>, Year.1 <dbl>, Year.2 <dbl>,
## #   duration <dbl>, LTER_site <chr>
ggplot(wrtds_trends, aes(estC))+
  geom_histogram()

ggplot(wrtds_trends, aes(lowC95))+
  geom_histogram()

ggplot(wrtds_trends, aes(upC95))+
  geom_histogram()

wrtds_trends %>% 
  filter(lowC95 < -2e7)
## # A tibble: 1 x 33
##   LTER  site  rejectC pValC   estC    lowC   upC lowC50  upC50  lowC95 upC95
##   <chr> <chr> <lgl>   <dbl>  <dbl>   <dbl> <dbl>  <dbl>  <dbl>   <dbl> <dbl>
## 1 UMR   M078~ FALSE   0.146 -0.446 -9.46e7 0.236 -0.871 -0.449 -9.53e7 0.906
## # ... with 22 more variables: likeCUp <dbl>, likeCDown <dbl>, rejectF <lgl>,
## #   pValF <dbl>, estF <dbl>, lowF <dbl>, upF <dbl>, lowF50 <dbl>, upF50 <dbl>,
## #   lowF95 <dbl>, upF95 <dbl>, likeFUp <dbl>, likeFDown <dbl>, baseConc <dbl>,
## #   baseFlux <dbl>, iBoot <dbl>, startSeed <dbl>, nBootGood <dbl>,
## #   Year.1 <dbl>, Year.2 <dbl>, duration <dbl>, LTER_site <chr>
wrtds_trends %>% 
  filter(upC95 > 2e7)
## # A tibble: 1 x 33
##   LTER  site  rejectC pValC  estC  lowC   upC lowC50 upC50 lowC95  upC95 likeCUp
##   <chr> <chr> <lgl>   <dbl> <dbl> <dbl> <dbl>  <dbl> <dbl>  <dbl>  <dbl>   <dbl>
## 1 ARC   Imna~ FALSE   0.795 0.553 -6.63  13.1 -0.595 0.718  -36.9 5.17e7      NA
## # ... with 21 more variables: likeCDown <dbl>, rejectF <lgl>, pValF <dbl>,
## #   estF <dbl>, lowF <dbl>, upF <dbl>, lowF50 <dbl>, upF50 <dbl>, lowF95 <dbl>,
## #   upF95 <dbl>, likeFUp <dbl>, likeFDown <dbl>, baseConc <dbl>,
## #   baseFlux <dbl>, iBoot <dbl>, startSeed <dbl>, nBootGood <dbl>,
## #   Year.1 <dbl>, Year.2 <dbl>, duration <dbl>, LTER_site <chr>
wrtds_trends_trimmed <- wrtds_trends %>% 
  filter(lowC95>-50) %>% 
  filter(upC95<50)
summary(wrtds_trends_trimmed)
##      LTER               site            rejectC            pValC        
##  Length:44          Length:44          Mode :logical   Min.   :0.01980  
##  Class :character   Class :character   FALSE:33        1st Qu.:0.09864  
##  Mode  :character   Mode  :character   TRUE :11        Median :0.23482  
##                                                        Mean   :0.33119  
##                                                        3rd Qu.:0.51405  
##                                                        Max.   :0.97723  
##                                                                         
##       estC              lowC              upC               lowC50        
##  Min.   :-1.5187   Min.   :-5.5108   Min.   :-0.01031   Min.   :-2.94079  
##  1st Qu.:-0.1363   1st Qu.:-0.5729   1st Qu.: 0.35689   1st Qu.:-0.25621  
##  Median : 0.1965   Median :-0.2876   Median : 0.86837   Median :-0.05223  
##  Mean   : 0.2363   Mean   :-0.6132   Mean   : 1.10630   Mean   :-0.08456  
##  3rd Qu.: 0.4931   3rd Qu.:-0.1114   3rd Qu.: 1.38779   3rd Qu.: 0.30061  
##  Max.   : 5.3402   Max.   : 0.9813   Max.   :10.16405   Max.   : 1.23391  
##                                                                           
##      upC50              lowC95             upC95             likeCUp       
##  Min.   :-0.58212   Min.   :-17.5133   Min.   : 0.02856   Min.   :0.02778  
##  1st Qu.:-0.01424   1st Qu.: -0.9011   1st Qu.: 0.43535   1st Qu.:0.25490  
##  Median : 0.42081   Median : -0.4742   Median : 1.13162   Median :0.71569  
##  Mean   : 0.55089   Mean   : -1.1159   Mean   : 1.37864   Mean   :0.60333  
##  3rd Qu.: 0.85097   3rd Qu.: -0.2064   3rd Qu.: 1.63045   3rd Qu.:0.94059  
##  Max.   : 6.97374   Max.   :  0.8973   Max.   :10.41049   Max.   :0.99505  
##                                                           NA's   :1        
##    likeCDown         rejectF            pValF             estF           
##  Min.   :0.004951   Mode :logical   Min.   :0.0198   Min.   :-10352.558  
##  1st Qu.:0.059406   FALSE:34        1st Qu.:0.1570   1st Qu.:    -0.056  
##  Median :0.284314   TRUE :10        Median :0.4122   Median :     0.033  
##  Mean   :0.396673                   Mean   :0.4393   Mean   :  -303.418  
##  3rd Qu.:0.745098                   3rd Qu.:0.6886   3rd Qu.:     6.917  
##  Max.   :0.972222                   Max.   :0.9244   Max.   :   278.731  
##  NA's   :1                                                               
##       lowF                 upF               lowF50          
##  Min.   :-20645.685   Min.   :   0.002   Min.   :-13299.324  
##  1st Qu.:   -16.438   1st Qu.:   0.133   1st Qu.:    -0.797  
##  Median :    -0.205   Median :   0.769   Median :    -0.002  
##  Mean   :  -697.895   Mean   : 302.531   Mean   :  -432.513  
##  3rd Qu.:    -0.002   3rd Qu.:  53.948   3rd Qu.:     0.241  
##  Max.   :    33.868   Max.   :5488.965   Max.   :   196.891  
##                                                              
##      upF50               lowF95               upF95              likeFUp       
##  Min.   :-4964.677   Min.   :-24081.227   Min.   :    0.002   Min.   :0.09406  
##  1st Qu.:    0.014   1st Qu.:   -37.228   1st Qu.:    0.176   1st Qu.:0.34259  
##  Median :    0.305   Median :    -0.325   Median :    0.915   Median :0.69444  
##  Mean   : -111.685   Mean   :  -816.794   Mean   :  738.278   Mean   :0.62305  
##  3rd Qu.:   17.943   3rd Qu.:    -0.012   3rd Qu.:   73.833   3rd Qu.:0.92574  
##  Max.   :  402.210   Max.   :     4.452   Max.   :18713.037   Max.   :0.99505  
##                                                               NA's   :3        
##    likeFDown           baseConc         baseFlux            iBoot       
##  Min.   :0.004951   Min.   : 0.661   Min.   :    0.03   Min.   : 50.00  
##  1st Qu.:0.074257   1st Qu.: 2.534   1st Qu.:    0.50   1st Qu.: 50.00  
##  Median :0.305556   Median : 4.725   Median :    3.16   Median : 50.00  
##  Mean   :0.376951   Mean   : 5.004   Mean   : 1717.12   Mean   : 61.93  
##  3rd Qu.:0.657407   3rd Qu.: 5.935   3rd Qu.:  204.11   3rd Qu.: 59.50  
##  Max.   :0.905941   Max.   :20.101   Max.   :42361.98   Max.   :100.00  
##  NA's   :3                                                              
##    startSeed        nBootGood          Year.1         Year.2    
##  Min.   :494817   Min.   : 50.00   Min.   :1964   Min.   :2012  
##  1st Qu.:494817   1st Qu.: 50.00   1st Qu.:1984   1st Qu.:2017  
##  Median :494817   Median : 50.00   Median :1995   Median :2018  
##  Mean   :494817   Mean   : 61.93   Mean   :1991   Mean   :2018  
##  3rd Qu.:494817   3rd Qu.: 59.50   3rd Qu.:1997   3rd Qu.:2019  
##  Max.   :494817   Max.   :100.00   Max.   :2003   Max.   :2019  
##                                                                 
##     duration      LTER_site        
##  Min.   :15.00   Length:44         
##  1st Qu.:21.75   Class :character  
##  Median :22.00   Mode  :character  
##  Mean   :27.00                     
##  3rd Qu.:31.50                     
##  Max.   :53.00                     
## 
ggplot(wrtds_trends_trimmed, aes(site, estC)) +
  geom_point()+
  geom_linerange(aes(ymin = lowC95, ymax = upC95))+
  coord_flip()+
  geom_hline(yintercept = 0)+
  scale_x_discrete(limits = rev)+
  facet_wrap(~LTER, scales='free')

#trying to make summary points per LTER didn't work this way

# ggplot(wrtds_trends_trimmed) +
#   geom_point(aes(site, estC))+
#   geom_linerange(aes(site, ymin = lowC95, ymax = upC95))+
#   coord_flip()+
#   geom_hline(yintercept = 0)+
#   scale_x_discrete(limits = rev)+
#   stat_summary(aes(LTER, fun = mean, geom = "point")) + 
#   stat_summary(aes(LTER, fun.data = mean_cl_normal,
#                geom = "pointrange", fun.args = list(mult = 1)))+
#   facet_wrap(~LTER, scales='free')

3.1 new attempt at weighted forest plots

# data stright from wrtds output

wrtds_trends_trimmed%>%
  ggplot(aes(fill = LTER, x = paste(LTER, site, sep = "::"), y = estC, ymin = lowC95, ymax=upC95))+
  geom_pointrange(color="black", shape = 21, size = 2)+
  scale_x_discrete(limits = rev)+
  coord_flip()+
  geom_hline(yintercept = 0, linetype="dashed")+
  theme(axis.text=element_text(size=16))+
  # scale_color_manual(values = jever2color)+
  # scale_fill_manual(values = jever2color)+
  #scale_size_manual(values=c(1,1.5))+
  #scale_shape_manual(values = c(21, 22))+
  xlab("site")+
  ylab("Slope of annual seasonality change")

#make SE data
wrtds_trends_trimmed <- wrtds_trends_trimmed %>% 
  mutate(SE_est = ((upC95-lowC95)/2)/sqrt(duration))

# do meta-analysis by LTER
mafits_trends <-  wrtds_trends_trimmed %>%
  group_by(LTER) %>% #group by type of nutrient manipulation and broad response
  do(model = tidy(rma(estC, SE_est, data = .), conf.int=TRUE)) %>% #do weighted meta-analysis for each grouping
  unnest(model)


#plot number observations and weighted data from above

num_obs <- wrtds_trends_trimmed %>% 
  group_by(LTER) %>%
  summarize(num_obs=length(estC), mean_estimate=mean(estC))
mafits_trends
## # A tibble: 9 x 9
##   LTER     term   type   estimate std.error statistic p.value conf.low conf.high
##   <chr>    <chr>  <chr>     <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
## 1 AND      overa~ summa~   0.0660    0.104      0.634 0.526     -0.138    0.270 
## 2 HBR      overa~ summa~   0.484     0.187      2.59  0.00973    0.117    0.850 
## 3 KRR      overa~ summa~   0.165     0.249      0.660 0.509     -0.324    0.654 
## 4 LMP      overa~ summa~   0.851     0.512      1.66  0.0964    -0.152    1.85  
## 5 LUQ      overa~ summa~   0.708     0.518      1.37  0.172     -0.308    1.72  
## 6 MCM      overa~ summa~  -0.113     0.0971    -1.16  0.244     -0.303    0.0773
## 7 NWT      overa~ summa~   2.50      2.52       0.991 0.322     -2.44     7.44  
## 8 UMR      overa~ summa~  -0.0829    0.357     -0.232 0.816     -0.783    0.617 
## 9 UMR_trib overa~ summa~   0.166     0.265      0.626 0.531     -0.354    0.686
num_obs
## # A tibble: 9 x 3
##   LTER     num_obs mean_estimate
##   <chr>      <int>         <dbl>
## 1 AND            7        0.0419
## 2 HBR            9        0.477 
## 3 KRR            4        0.167 
## 4 LMP            1        0.851 
## 5 LUQ            2        0.372 
## 6 MCM            8       -0.122 
## 7 NWT            2        2.80  
## 8 UMR            5       -0.164 
## 9 UMR_trib       6       -0.0415
mafits_trends <- left_join(mafits_trends, num_obs)

mafits_trends%>%
  ggplot(aes(x = LTER, y = estimate, ymin = conf.low, ymax=conf.high))+
  geom_pointrange(color="black", shape = 21, size = 1.8, fill = "#218D3A")+
  scale_x_discrete(limits = rev)+
  coord_flip()+
  geom_hline(yintercept = 0, linetype="dashed")+
  geom_text(size=4, fontface="italic", position = position_nudge(x = 0.35, y=.25),
            aes(label = num_obs))+
  theme(axis.text=element_text(size=16))+
  scale_color_manual(values = jever2color)+
  scale_fill_manual(values = jever2color)+
  #facet_wrap(~variable)+
  #scale_size_manual(values=c(1,1.5))+
  #scale_shape_manual(values = c(21, 22))+
  xlab("LTER")+
  ylab("Slope DSi change over time")

#ggsave("wrtds_trends_weighted.png")

3.2 old attempts, maybe not useful

calculating a rough estimate of standard error from the 95% CIs for now, but need to make sure this is ok to do (should still give a relatively good weighting for meta-regression stuff for now/makes sense when compared to the CIs)

wrtds_trends_trimmed <- wrtds_trends_trimmed %>% 
  mutate(SE_est = ((upC95-lowC95)/2)/sqrt(duration))

ggplot(wrtds_trends_trimmed, aes(site, SE_est))+
  geom_point()+
  coord_flip()

trying to do a meta-regression model for each LTER because that’s the only way I know how to get overall weighted effect sizes for each LTER, but metafor doesn’t like whatever I’m doing

likelihood vs pval stuff to check out how those are all related

wrtds_trends %>% count(rejectC)
## # A tibble: 2 x 2
##   rejectC     n
##   <lgl>   <int>
## 1 FALSE      39
## 2 TRUE       13
ggplot(wrtds_trends, aes(pValC, rejectC))+
  geom_point()+
  geom_vline(xintercept=.05)

ggplot(wrtds_trends, aes(likeCDown, rejectC))+
  geom_point()+
  geom_vline(xintercept=.05)

ggplot(wrtds_trends, aes(likeCUp, rejectC))+
  geom_point()+
  geom_vline(xintercept=.95)

4 maps and latitude (non)effects

importing lat-long data and having a quick look- looks like Sagehen needs to swap hemispheres

latlongs <- read_csv("LongTermWatersheds_LatLong.csv") %>% 
  rename(site=Stream.Site)

leaflet(latlongs) %>%
  addTiles() %>% 
  addMarkers(lng = ~Longitude, lat = ~Latitude, popup = ~ LTER)
latlongs <- latlongs %>% 
  mutate(Longitude = ifelse(LTER=='Sagehen', -Longitude, Longitude))

leaflet(latlongs) %>%
  addTiles() %>% 
  addMarkers(lng = ~Longitude, lat = ~Latitude, popup = ~ LTER)
a <- unique(wrtds_trends$site) #usi this to check where there are differences in the colnames
b <- unique(latlongs$site)
setdiff(b,a)
## [1] "TW Weir"      "Toolik Inlet" "Q1"           "Q2"           "Q3"          
## [6] "QP"
b %in% a
##  [1]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE FALSE  TRUE  TRUE
## [13]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [25]  TRUE  TRUE  TRUE FALSE FALSE FALSE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE
## [37]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [49]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
which(a %in% b)
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
## [51] 51 52
wrtds_trends2 <- left_join(wrtds_trends, latlongs, by='site')
wrtds_slopes_ES_Si2 <- left_join(wrtds_slopes_ES_Si, latlongs, by="site")
wrtds_slopes_ES_discharge2 <- left_join(wrtds_slopes_ES_discharge, latlongs, by="site")
wrtds_slopes_ES_flux2 <- left_join(wrtds_slopes_ES_flux, latlongs, by="site")

joining with slope data- here these represent untransformed values of the slopes (“real” slope units of change per time from WRTDS trend estimates I think) where all above forest plots were Z-transformed

pal <- colorBin(palette = "RdBu", bins = c(-6,-3,-2,-1,0,1,2,3,6), domain = wrtds_trends$estC, reverse = TRUE)

#pal <- colorNumeric("RdBu", wrtds_trends$estC)

leaflet(wrtds_trends2) %>%
  addTiles() %>% 
  addCircleMarkers(lng = ~Longitude,
             lat = ~Latitude,
             #popup = ~ LTER.x,
             color="black",
             stroke=TRUE,
             fillColor = ~pal(estC),
             fillOpacity = 10,
             opacity = .5,
             label = ~paste0(wrtds_trends2$LTER.x,"::", wrtds_trends2$site, ", slope: ", estC),
             radius = 10) 
# I don't think the legend is super helpful..

  # addLegend(position = "bottomleft",
  #                 pal = pal,
  #                 values = range(wrtds_trends$estC),
  #                 title = "WRTDS slope estimate (DSi~year)")

4.1 meta-regression including absolute latitude as variable

plotting absolute lat vs slopes for DSi/discharge/flux

ggplot(wrtds_trends2, aes(abs(Latitude), estC))+
  geom_point()+
  geom_smooth(method='lm')

ggplot(wrtds_slopes_ES_Si2, aes(abs(Latitude), estimate))+
  geom_point()+
  geom_smooth(method='lm')+
  ylab('slope DSi from annual estimates')

ggplot(wrtds_slopes_ES_discharge2, aes(abs(Latitude), estimate))+
  geom_point()+
  geom_smooth(method='lm')+
  ylab('slope discharge from annual estimates')

ggplot(wrtds_slopes_ES_flux2, aes(abs(Latitude), estimate))+
  geom_point()+
  geom_smooth(method='lm')+
  ylab('slope DSi flux from annual estimates')

meta-regression stuff also shows no effect of lat, but here’s where we could also plug in additional variables to see if anything modifies slopes

mod_Si3<-rma.mv(yi,vi,
          data=wrtds_slopes_ES_Si2,
          random = ~1|site,
          mods = ~abs(Latitude))
summary(mod_Si3)
## 
## Multivariate Meta-Analysis Model (k = 52; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
## -47.0893   94.1786  100.1786  105.9147  100.7004   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed  factor 
## sigma^2    0.3424  0.5851     52     no    site 
## 
## Test for Residual Heterogeneity:
## QE(df = 50) = 512.2841, p-val < .0001
## 
## Test of Moderators (coefficient 2):
## QM(df = 1) = 3.8262, p-val = 0.0505
## 
## Model Results:
## 
##                estimate      se     zval    pval    ci.lb   ci.ub 
## intrcpt          0.5588  0.2682   2.0834  0.0372   0.0331  1.0846  * 
## abs(Latitude)   -0.0108  0.0055  -1.9561  0.0505  -0.0217  0.0000  . 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod_discharge2<-rma.mv(yi,vi,
          data=wrtds_slopes_ES_discharge2,
          random = ~1|site,
          mods = ~abs(Latitude))
summary(mod_discharge2)
## 
## Multivariate Meta-Analysis Model (k = 52; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
##   1.5312   -3.0624    2.9376    8.6737    3.4593   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed  factor 
## sigma^2    0.0105  0.1023     52     no    site 
## 
## Test for Residual Heterogeneity:
## QE(df = 50) = 64.3719, p-val = 0.0832
## 
## Test of Moderators (coefficient 2):
## QM(df = 1) = 1.4759, p-val = 0.2244
## 
## Model Results:
## 
##                estimate      se     zval    pval    ci.lb   ci.ub 
## intrcpt          0.2798  0.1014   2.7599  0.0058   0.0811  0.4785  ** 
## abs(Latitude)   -0.0026  0.0021  -1.2149  0.2244  -0.0067  0.0016     
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod_flux2<-rma.mv(yi,vi,
          data=wrtds_slopes_ES_flux2,
          random = ~1|site,
          mods = ~abs(Latitude))
summary(mod_flux2)
## 
## Multivariate Meta-Analysis Model (k = 52; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
##  -4.9865    9.9729   15.9729   21.7090   16.4947   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed  factor 
## sigma^2    0.0316  0.1777     52     no    site 
## 
## Test for Residual Heterogeneity:
## QE(df = 50) = 88.8180, p-val = 0.0006
## 
## Test of Moderators (coefficient 2):
## QM(df = 1) = 2.5593, p-val = 0.1096
## 
## Model Results:
## 
##                estimate      se     zval    pval    ci.lb   ci.ub 
## intrcpt          0.3818  0.1195   3.1959  0.0014   0.1476  0.6159  ** 
## abs(Latitude)   -0.0040  0.0025  -1.5998  0.1096  -0.0088  0.0009     
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1